function SpectralClustNG()
%% documentation
% takes as input data in the form of a matrix
% a value for sigma
% a value for k (the number of clusters)

%% dataset1
cluster1 = [2.49671415,3.57921282;1.8617357,2.76743473;2.64768854,1.53052561;...   
            3.52302986,2.54256004;1.76584663,1.53658231;1.76586304,1.53427025];    
cluster2 = [ 7.24196227,6.09197592;5.08671976,5.5876963;5.27508217,8.46564877;...
           6.43771247,6.7742237;5.98716888,7.0675282;7.31424733,5.57525181];    
data1 = vertcat(cluster1,cluster2);
sizeC1 = length(cluster1);
sizeC2 = length(cluster2);
labels1 = vertcat(zeros(sizeC1,1)+1,zeros(sizeC2,1)+2);

%% dataset2
cluster1 = [1.10,3.3;1.18,3.3;1.26,3.3;1.26,3.2;1.22,3.1;1.14,3.1;1.06,3.1;0.98,3.1;0.90,3.1;0.82,3.1;0.75,3.15;0.71,3.2;...
            0.7,3.3;0.7,3.4;0.7,3.5;0.70,3.6;0.70,3.7;0.73,3.8;0.8,3.88;0.9,3.9;1.0,3.9;1.07,3.88;1.15,3.83;1.21,3.75;1.26,3.65];

cluster2 = [1.7,3.9;1.7,3.8;1.7,3.7;1.7,3.6;1.7,3.5;1.7,3.4;1.7,3.3;1.7,3.2;1.7,3.1;1.76,3.8;1.82,3.7;1.88,3.6;1.94,3.5;...
            2.00,3.4;2.06,3.3;2.12,3.2;2.18,3.1;2.25,3.9;2.25,3.8;2.25,3.7;2.25,3.6;2.25,3.5;2.25,3.4;2.25,3.3;2.25,3.2;2.25,3.1];

cluster3 = [2.7,3.9;2.7,3.8;2.7,3.7;2.7,3.6;2.7,3.5;2.7,3.4;2.7,3.3;2.7,3.2;2.74,3.15;2.80,3.1;2.90,3.1;2.98,3.1;3.06,3.1;3.16,3.1;...                     
           3.25,3.9;3.25,3.8;3.25,3.7;3.25,3.6;3.25,3.5;3.25,3.4;3.25,3.3;3.25,3.2;3.21,3.15];  
        
data2 = vertcat(cluster1,cluster2,cluster3);
sizeC1 = length(cluster1);
sizeC2 = length(cluster2);
sizeC3 = length(cluster3);
labels2 = vertcat(zeros(sizeC1,1)+1,zeros(sizeC2,1)+2,zeros(sizeC3,1)+3);

%% specify which data to use and k
k = 3;
sigma = 2.0;
mat = data2;
labels = labels2;
[n,d] = size(mat);


%% create the similariyt matrix from raw matrix
values = pdist(mat,'euclidean');
sMat = squareform(values);
size(sMat)

%% create affinity matrix from distance matrix
aMat = exp(- sMat.^2 / 2 * sigma^2);

%% create the diagonal matrix
D = diag(sum(aMat, 2))^ (-0.5);

%% create the L matrix
L = D * aMat * D;
size(L)

%% create the x matrix
[eigVecs, eigVals] = eigs(L,k);
xMat = eigVecs;

%% create the normalized y matrix
unitLengths = sum(xMat.^2,1).^(0.5);
yMat = zeros(n,k);

for col = (1:k)
    yMat(:,col) = xMat(:,col) / unitLengths(col);    
end

%% here we use the result from eigs or above
[IDX,C] = kmeans(yMat,k,'replicates',5);
horzcat(IDX,labels);

%% make plot
figure()
hold on
x = mat(:,1);
y = mat(:,2);
colors = ['b','k','r','g'];
for i = 1:k
    scatter([x(IDX==i)],[y(IDX==i)],'ks','MarkerFaceColor',colors(i))
end
buffer = 0.5;
xlim([min(x) - buffer, max(x) + buffer]);
ylim([min(y) - buffer, max(y) + buffer]);
hold off

%% debugging
%sMat(1:10,1:10)
%aMat(1:10,1:10)
%L(1:10,1:10)
xMat(1:10,1:k)
yMat(1:10,1:k)

[eVecs,eVals] = eig(L);
eVecs = eVecs * -1.0;
eVecs(1:5,:)
sum(abs(eVecs),1)


